NASA 

Technical 

Paper 

214L_ 

sr- . T ■’ir. •. 

April 1992 


•_4 .. : .. > . , >=.- *~ T - ' ‘ 




r 



An Efficient HZETRN 
(A Galactic Cosmic 
Ray Transport Code) 


Judy L, Shinn 
and John W. Wilson 



(\A$A-TP-3147) AN Iff ICIENT HZET^n (A 
GALACTIC COSMIC KAY TRANSPORT CODE) (NASA) 
17 p CSCL 03B 



Unci as 

Hl/93 0084320 



NASA 

Technical 

Paper 

3147 


1992 


An Efficient HZETRN 
(A Galactic Cosmic 
Ray Transport Code) 


Judy L. Shinn 
and John W. Wilson 
Langley Research Center 
Hampton , Virginia 


NASA 

National Aeronautics and 
Space Administration 

Office of Management 

Scientific and Technical 
Information Program 



Abstract 


To meet the challenge of the future deep-space program, which involves 
extended manned space missions, an accurate and efficient engineering code 
for analyzing the shielding requirement against the high-energy galactic heavy 
ions is needed. The HZETRN is a deterministic code developed at Langley 
Research Center that is constantly under improvement both in physics and 
numerical computation and is targeted for such use. One problem area con- 
nected with the space-marching technique used in this code is the propagation 
of the local truncation error . By improving the numerical algorithms for in- 
terpolation, integration, and grid distribution formula, the efficiency of the 
code is increased by a factor of eight as the number of energy grid points is re- 
duced. The numerical accuracy of better than 2 percent for a shield thickness 
of 150 g/cm? is found when a 45-point energy grid is used. The propagating 
step size, which is related to the perturbation theory, is also reevaluated. 


Introduction 

As the space program proceeds into an era of 
extended manned space operations, the shielding 
from galactic heavy ions becomes a problem of ever- 
increasing importance (ref. 1). The high-energy 
heavy ions originating in deep space interact with 
target nuclei resulting in energy degradation and nu- 
clear fragmentations. These fragmentations produce 
secondary and subsequent-generation reaction prod- 
ucts that alter the elemental and isotopic composi- 
tion of the transported radiation fields. A realistic 
estimate of flux in a critical organ of interest can be 
made when the nuclear fragmentation data become 
available as inputs to the galactic cosmic ray trans- 
port code (HZETRN (ref. 2)), recently developed at 
Langley Research Center. 

As NASA places efforts on the experimental pro- 
gram to produce fragmentation data, the space ra- 
diation transport codes including HZETRN are fur- 
ther being improved, refined, and updated to meet 
future mission requirements. One such effort is to im- 
prove the efficiency and numerical accuracy of these 
deterministic codes, such that the codes can easily 
be used as engineering tools and will also accom- 
modate future expansion for more sophistication in 
physics. Recently, the work on the baryon trans- 
port code (BRYNTRN) is an example in this area 
(ref. 3). This code, as well as HZETRN, is based on 
a space-marching formalism that calls special atten- 
tion to the error propagation. An analysis made in 
the study (ref. 3) showed that the propagated error 
tends to grow with the marching steps to a maxi- 
mum, but is proportional to the local error that can 


be minimized by improving various numerical algo- 
rithms and the grid generation. The results of these 
modifications have substantially improved the accu- 
racy and efficiency of BRYNTRN. 

Because HZETRN calculates the transport of the 
galactic heavy ions through target materials, the 
number of species and the range of energy considered 
in this code are much larger than in BRYNTRN. For 
this reason, the benefit from the improvement of nu- 
merical computation in HZETRN is expected to be 
far greater than in BRYNTRN. In this report, the 
description and testing of the changes to the compu- 
tational procedures in HZETRN as well as heavy-ion 
transport theory are presented. The improvements 
in efficiency and accuracy are also evaluated. 

Galactic Cosmic Ray Transport Method 

Galactic Cosmic Ray Transport Theory 

In moving through extended matter, heavy ions 
lose energy through interaction with atomic electrons 
along their trajectories. On occasion, they interact 
violently with nuclei of the matter, producing ion 
fragments moving in the forward direction and low- 
energy fragments of the struck target nucleus. The 
transport equations for the short-range target frag- 
ments can be solved in closed form in terms of colli- 
sion density (ref. 4). Hence, the projectile fragment 
transport is the interesting unsolved problem. In pre- 
vious work, the projectile ion fragments were treated 
as if all went straight ahead (ref. 5). The straight- 
ahead approximation is found to be quite accurate 
for the nearly isotropic cosmic ray fluence (ref. 4). 



With the straight-ahead approximation and the target secondary fragments neglected (refs. 4 and 5), the 
transport equation may be written as 
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where 3>j(x, E ) is the flux of ions of type j with atomic mass Aj at x moving along the x-axis at energy E in 
units of MeV/amu, <jj is the corresponding macroscopic nuclear absorption cross section, Sj(E) is the change 
in E per unit distance, and mj ^ is the multiplicity of ion j produced in collision by ion k. The corresponding 
nucleon transport equation (refs. 3, 6, and 7) is 
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The rrijk , aj are assumed energy independent in equation (1) and the full energy dependence is retained in 
equation (2). The solution to equations (1) and (2) is to be found subject to boundary specification at x = 0 
and arbitrary E as 

,E) = Fj{E) (3) 

Usually Fj(E) is called the incident beam spectrum. 

It follows from Bethe’s theory (ref. 8) that 
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and holds for all energies above 100 keV/amu provided the ions remain fully stripped. The range of the ion is 
given as 
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The subscript p refers to proton. Equation (6) is quite accurate at high energy and only approximately true at 
low energy because of electron capture by the ion that effectively reduces its charge (ref. 9), higher order Born 
corrections to Bethe’s theory (ref. 10), and nuclear stopping at the lowest energies (refs. 11 and 12). Herein, 
the parameter i/j is defined as 

VjRj{E) = v k R k {E) (7) 


so that 
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Equations (6) to (8) are used in the subsequent development and the energy variation in vj is neglected. 
A method of solution is now discussed. For the purpose of solving equation (1), define the coordinates 
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where rjj varies along the particle path and is constant along the particle trajectory. The new fluence 
functions are taken as 


Xj{Vj-Aj) = Sj{E)(f>j{x,E) - rj) 
XkiVj^j) = Xk(Vk^k) 


where 
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and rj = Rj(E). Under this coordinate mapping, equation (1) becomes 
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where crj are assumed to be energy 7 independent. There is a small variation in (^20 percent), which must 
eventually be taken into account. Solving equation (15) by using line integration with an integrating factor 


results in 


dj (Vj ,Zj) = exp - tTj (£, + V ] ) 


Xj(Vj,Zj) = exp - - + rjj) Xj{-Zj,Zj) 


(16) 


1 Hi 

+ 2./, eXP 


\ °jiv' -Vj) m J k a k 77 XfcWUfc) 

^ J k Vk 


(17) 


where 


, v k + 7 , Vk - Vj, , ,/ v k - Vj v k + V 

* = — + and = -2^ + 


2^Jt 


With equation (11) one may show 
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where h is the step size in the x direction. 

It is clear from equation (18) that 
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which is correct to order h 2 . This expression may be further approximated by 
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which is accurate to — vj ) /i j . Equation (22) is the basis of the galactic cosmic ray (GCR) transport 

code GCRTRN (refs. 13 to 15). A few years ago, GCRTRN and the baryon transport code (BRYNTRN) were 
coupled together as a new code (HZETRN) which effectively solves equation (2) by adding a heavy ion collision 
source of nucleons to the right-hand side of the equation. Equation (22) provides the propagating algorithm 
for the heavy ions. The corresponding propagating procedure for the nucleons is given as (refs. 3 and 6) 
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with the order of h 2 . 

There are several quantities of interest that are now given. The integral fluence is given as 
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The energy absorption per gram is 
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with the dose equivalent given as 
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where Qp is the quality factor. These quantities are used in shield design studies for protection against galactic 
cosmic rays. 


Numerical Procedure 

The secondary particle production term of the propagation algorithm for nucleons in equation (23) has been 
further reduced to a form that can be implemented with ease for numerical integration. The details of the 
form and its validity have been discussed elsewhere (ref. 6) and will not be repeated here. For the heavy ions, 
the secondary production term (the second term on the right side of eq. (22)) does not involve any integration, 
however, the interpolation of the transformed fluence function is based on the independent variable r fc , which is 
different from ry, the range of ions of type j given at the left side of the equation. To circumvent the problem, 
the equation is further modified. Recall the definition of Sj(E) with E given in units of MeV/amu 
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with 

Sj(Ej) = Z] SpiEj/Aj) (28) 

where S p is the proton stopping power and Ej is the energy in MeV of ions of type j. It follows that 

Sj(E) = - -Z] S p (Ej/Aj) = VjSp(Ep) (29) 
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where 

Rewriting equation (11) as 
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we can define the new fluence function 
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Equation (22) now becomes 


E dE 1 


=/v 

Jo bp 


(E') 


ipj(x + h,r) — e a i h ipj{x, r + vjh) + ^ m jk a k 


o h g 


a k a j / u j 


— i!)' k {x,r + vjh) 


(30) 

(31) 

(32) 

(33) 

(34) 


so that there is only one single definition of range that is related to energy. The equation can now be solved 
by setting up the r (proton range) grid and marching the solution from x — 0 by steps of h to the desired 
thickness. 


Error Propagation 

In considering how errors are propagated in the use of equation (34), the error is introduced locally by 
calculating ^(x,r + Vjh) over the range (energy) grid. Limiting our current analysis to the first term of 
equation (34), it is defined at each range grid r* that 

H- h,rj) ~ e~°i h ^'(x,r* + Vjh) (35) 

We denote the truncation error e l introduced in the interpolation procedure to the interpolated value, 
as 
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because hoj 1. Clearly the propagated error on the mth step is bound by 
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where e(h) is the maximum error per step. With the 
increasing value of ra, the propagated error grows 
each step to a maximum value of e(h)/haj. Because 
the increase of h value is limited by the perturbation 
theory, reducing the local truncation error is the only 
viable approach left for reducing the propagated error 
to a desired level. The same consideration may be 
applied for the second term of equation (34), as the 
terms are of similar nature. 

Numerical Algorithms 

The error analysis shown in the previous section 
has concluded that to effectively reduce the level of 
propagated error, the local truncated error must be 
reduced. There are three basic numerical algorithms 
that are involved in solving equations (23) to (2G) and 
(34): interpolation, integration, and grid generation. 
The integration scheme does not affect the error 
propagation for the heavy-ion transport, but does 
affect the nucleons and the dose calculations. The 
choice of a grid distribution that is interrelated to 
the interpolation and integration scheme can increase 
the efficiency of the code if the number of grids can 
be reduced. 

The interpolation scheme to be used here is the 
third-order Lagrange method as used successfully in 
the work for improving BRYNTRN (ref. 3). With 
the four neighboring interpolating grids (data points) 
placed evenly on both sides of the interpolated point, 
the error will tend to be the smallest in the middle 
interval of all the data points if the grid distribution 
is rather uniform (ref. 16). The choice of a much 
higher order Lagrange method will substantially de- 
crease the efficiency of the code, because there are 
more than 10 interpolation calls for each single en- 
ergy point at each step. Other interpolation methods 
such as a cubic spline were considered but discarded. 
The splines are, in general, more accurate. However, 
their characteristic large excursions (oscillations) can 
result in erroneous, unpredictable solutions. 

The same procedure for numerical integration 
used for the improved BRYNTRN (ref. 3) will also 
be used here for HZETRN. It is based on the com- 
pound quadrature formulation summing over all the 
subintervals between the grids with the midpoint 
evaluated by making use of the improved interpola- 
tion procedure mentioned above. A simple numerical 
method such as Simpson’s rule is used to integrate for 
the subintervals. 

There arc three binding conditions that dictate 
how the grids should be distributed. The first is the 
shape of the input spectrum. Because the galactic 
cosmic ray fluences are several orders of magnitude 


larger at the lower energy end (ref. 17), the loga- 
rithmic scale will be used for the energy or range 
coordinate as was done in reference 3. The second 
condition is related to the choice of the interpolation 
method that requires the four neighboring grids to 
be as uniformly spaced (on logarithm scale) as possi- 
ble so that the interpolation error can be minimized. 
Because the interpolation is performed on the range 
grid rather than on the energy grid, a uniform grid 
distribution on a logarithm of range r is desired. The 
third is related to the efficiency of the code that is 
found to increase almost quadratically with the de-' 
crease of grid points. With the uniform grid points 
on logarithm r scale as the basic structure, the distri- 
bution can further be modified to reduce the number 
of points in the region in which the information is not 
propagating through the steps. For BRYNTRN, it is 
the region below r m * n + /i, or approximately 1 g/cm 2 
because r m \ n 1 and assuming h — 1 g/cm 2 (ref. 3). 
The same applies for HZETRN, although the inter- 
polation is now at r m j n -f- Vjh, where uj is always 
equal to or much greater than 1. 

Results and Discussion 

In a previous study (ref. 3), the new grid distribu- 
tion with the number of grid points N equals 30 was 
found satisfactory for BRYNTRN, with the calcu- 
lated dose of 5 percent accuracy for a shield thickness 
of 150 g/cm 2 . Because the galactic cosmic rays are 
much harder than the solar flare protons, the upper 
limit of the energy range is usually taken to be about 
50 GeV as compared with a few GeV for solar pro- 
tons. Thus, more grid points are used for HZETRN. 
Tests were performed to determine the differences be- 
tween the interpolated ^(0, r + Vjh) and the analyt- 
ical results using the new interpolation method and 
grid distribution with N — 45, for all the f s and fc’s 
wdiere j < k. Samples of the results are displayed in 
figures 1 to 6. The overall error has been found to be 
less than 0.2 percent, with the particular grid gener- 
ation formula purposely adjusted so that the larger 
errors are absorbed at the low-energy region for r less 
than r m j n T v } h, where h = 1 g/cm 2 . The errors re- 
ported here are far less than the interpolation error 
reported in reference 3 for BRYNTRN. 

To test the convergence of the solution from the 
improved HZETRN (including the new integration 
procedure), the absorbed doses in tissue behind var- 
ious thicknesses of aluminum shield exposed to the 
galactic cosmic rays at solar minimum are calculated 
with N = 45, 60, and 90. The shield thickness 
is varied up to 150 g/cm 2 as indicated in figure 7. 
The doses from the heavy ions decrease significantly 
through the shield as they are fragmented by the 
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target nuclei. The doses from the nucleons increase 
rapidly to a maximum as a result of such fragmenta- 
tion and then gradually decrease. The relative errors 
in the absorbed dose for 45 grids or 60 grids compared 
with the dose for 90 grids are plotted in figures 8(a) 
and 8(b). The error for the proton dose increases 
with the shield thickness (see fig. 8(a)) as was ex- 
pected from the analysis of error propagation. There 
is some indication of oscillation for the neutron com- 
ponent (see fig. 8(b)), which was pointed out earlier 
(ref. 3) to be the result of the rapidly varying cross- 
section data for neutrons at low energy. The errors 
for the other component are found to be insignifi- 
cant compared with those for the nucleons. Thus, 
the maximum error is 2 percent for 45 grids and de- 
creases to less than 1 percent for 60 grids, showing 
good convergence. 

Another convergence issue that needs to be ad- 
dressed is the perturbation theory used in the trans- 
port theory. The perturbation theory requires 
( 7jh 1. Because uj is on the order of 0.01, we usu- 

ally take h = 1 g/cm 2 . Table 1 shows the effect 
of step size on some of the calculated doses, with 
h = 0.5 and 1 g/cm 2 . The differences between these 
two step sizes are insignificantly small, therefore, 
h — 1 g/cm 2 is retained for HZETRN. The overall 
efficiency for this code is improved about eight times 
as the grid points are reduced from 160 to 45, with an 
accuracy of 2 percent for the large shield thickness. 
A typical run time for producing the results shown 
in figure 7 is about 5000 sec on a CYBER 800 series 
computer. 

Concluding Remarks 

The efficiency of HZETRN (a galactic cosmic ray 
transport code) has been improved by approximately 
a factor of eight. The numerical algorithms for 
interpolation, integration, and energy grid generation 
were modified such that the number of grid points 
needed was reduced from 160 to 45 points. The error 
in dose calculation for 45 grid points was determined 
to be within 2 percent by a convergence test. The 
adequacy of the step size, which is related to the 
perturbation theory, was also examined. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
February 4, 1992 
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Table 1. Effect of Step Size on Calculated Doses Through 
Various Thicknesses of Aluminum Shield 


Aluminum 

shield 

thickness, 

g/cm 2 

Neutron dose, Gy 

Proton dose, Gy 

Alpha particle dose, Gy 

Step size 
= 0.5 g/cm 2 

Step size 
= 1 g/cm 2 

Step size 
= 0,5 g/cm 2 

Step size 
= 1 g/cm 2 

Step size 
= 0.5 g/cm 2 

Step size 
= 1 g/cm 2 

5 

6.3599E-3 

6.3315E-3 

8.2623E-2 

8.2375E-2 

2.4349E-2 

2.4349E-2 

10 

1.1232E-2 

1.1181E-2 

9.1011E-2 

9.0682E-2 

2.1400E-2 

2.1400E-2 

15 

1.4942E-2 

1.4873E-2 

9.4520E-2 

9.4148E-2 

1.8784E-2 

1.8784E-2 

20 

1.7738E-2 

1.7655E-2 

9.5463E-2 

9.5068E-2 

1.6499E-2 

1.6499E-2 




Figure 1. Fractional error of interpolating alpha-particle component of GCR spectrum at solar minimum with 
new grid distribution and interpolation method for generation of secondary alpha-particles. 



Figure 2. Fractional error of interpolating phosphorus component of GCR spectrum at solar minimum with 
new grid distribution and interpolation method for generation of secondary alpha-particles. 
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Figure 3. Fractional error of interpolating nickel component of GCR spectrum at solar minimum with new 
grid distribution and interpolation method for generation of secondary alpha-particles. 



Figure 4. Fractional error of interpolating phosphorus component of GCR spectrum at solar minimum with 
new grid distribution and interpolation method for generation of secondary phosphorus. 
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Figure 5. Fractional error in interpolating nickel component of GCR spectrum at solar minimum with new 
grid distribution and interpolation method for generation of secondary phosphorus. 



Figure 6. Fractional error in interpolating nickel component of GCR spectrum at solar minimum with new 
grid distribution and interpolation method for generation of secondary nickel. 
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Figure 7. Calculated absorbed dose components in tissue as function of aluminum shield thickness exposed to 
GCR at solar minimum. 
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Figure 8. Relative errors in doses compared with those calculated with 90 
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